function normal_equations()
% NORMAL_EQUATIONS 正规方程法求解最小二乘问题
% 
% 演示基于正规方程的最小二乘求解方法及其优缺点

fprintf('=== 正规方程法 ===\n\n');

%% 1. 正规方程的推导
fprintf('1. 正规方程的推导\n');
fprintf('最小二乘问题: min ||Ax - b||₂²\n');
fprintf('正规方程: A^T*A*x = A^T*b\n\n');

% 简单示例
A = [1, 1;
     1, 2;
     1, 3;
     1, 4];
b = [2.1; 2.9; 4.2; 4.8];

fprintf('系数矩阵 A (%d×%d):\n', size(A));
disp(A);
fprintf('观测向量 b:\n');
disp(b');

% 构造正规方程
AtA = A' * A;
Atb = A' * b;

fprintf('正规方程系数矩阵 A^T*A:\n');
disp(AtA);
fprintf('右端向量 A^T*b:\n');
disp(Atb');

% 求解正规方程
x_normal = AtA \ Atb;
fprintf('正规方程解: x = [%.6f; %.6f]\n', x_normal);

% 计算残差
residual = A * x_normal - b;
residual_norm = norm(residual);
fprintf('残差范数: ||r|| = %.6f\n', residual_norm);

%% 2. 几何解释
fprintf('\n2. 几何解释\n');
fprintf('正规方程的解使得残差向量与A的列空间正交\n\n');

% 验证正交性条件
orthogonality_check = A' * residual;
fprintf('正交性检验 A^T*r:\n');
disp(orthogonality_check');
fprintf('||A^T*r|| = %e (应该接近0)\n', norm(orthogonality_check));

% 投影解释
P = A * (AtA \ A');  % 投影矩阵
b_proj = P * b;      % b在列空间上的投影

fprintf('投影矩阵性质验证:\n');
fprintf('P² - P 的范数: %e\n', norm(P^2 - P, 'fro'));
fprintf('P^T - P 的范数: %e\n', norm(P' - P, 'fro'));

%% 3. 条件数问题
fprintf('\n3. 条件数问题\n');
fprintf('A^T*A的条件数是A的条件数的平方\n\n');

cond_A = cond(A);
cond_AtA = cond(AtA);

fprintf('cond(A) = %.2e\n', cond_A);
fprintf('cond(A^T*A) = %.2e\n', cond_AtA);
fprintf('cond(A)² = %.2e\n', cond_A^2);
fprintf('条件数平方关系的误差: %.2e\n', abs(cond_AtA - cond_A^2));

%% 4. 病态问题示例
fprintf('\n4. 病态问题示例\n');
fprintf('当A接近列相关时，A^T*A变得极度病态\n\n');

% 构造病态矩阵
epsilon = 1e-8;
A_ill = [1, 1;
         1, 1+epsilon;
         1, 1+2*epsilon;
         1, 1+3*epsilon];
b_ill = [1; 2; 3; 4];

fprintf('病态矩阵 A (ε = %e):\n', epsilon);
disp(A_ill);

AtA_ill = A_ill' * A_ill;
Atb_ill = A_ill' * b_ill;

cond_A_ill = cond(A_ill);
cond_AtA_ill = cond(AtA_ill);

fprintf('条件数比较:\n');
fprintf('cond(A) = %.2e\n', cond_A_ill);
fprintf('cond(A^T*A) = %.2e\n', cond_AtA_ill);

% 求解病态正规方程
x_ill = AtA_ill \ Atb_ill;
fprintf('病态正规方程解: x = [%.6f; %.6f]\n', x_ill);

% 与直接最小二乘比较
x_direct = A_ill \ b_ill;
fprintf('直接最小二乘解: x = [%.6f; %.6f]\n', x_direct);
fprintf('解的差异: %e\n', norm(x_ill - x_direct));

%% 5. 扰动分析
fprintf('\n5. 扰动分析\n');
fprintf('研究数据扰动对解的影响\n\n');

% 给观测数据添加小扰动
delta_b = 1e-6 * randn(size(b_ill));
b_perturbed = b_ill + delta_b;

% 扰动后的解
Atb_perturbed = A_ill' * b_perturbed;
x_perturbed = AtA_ill \ Atb_perturbed;

% 误差分析
rel_error_b = norm(delta_b) / norm(b_ill);
rel_error_x = norm(x_perturbed - x_ill) / norm(x_ill);

fprintf('数据相对扰动: %e\n', rel_error_b);
fprintf('解的相对误差: %e\n', rel_error_x);
fprintf('误差放大因子: %.2e\n', rel_error_x / rel_error_b);
fprintf('理论上界 κ(A^T*A): %.2e\n', cond_AtA_ill);

%% 6. Cholesky分解求解
fprintf('\n6. Cholesky分解求解正规方程\n');
fprintf('当A^T*A正定时，可以使用Cholesky分解\n\n');

% 使用良态矩阵
A_good = randn(10, 4);
b_good = randn(10, 1);
AtA_good = A_good' * A_good;
Atb_good = A_good' * b_good;

fprintf('良态矩阵条件数: %.2e\n', cond(AtA_good));

% 检查正定性
eigenvals = eig(AtA_good);
fprintf('A^T*A的最小特征值: %.6e\n', min(eigenvals));
fprintf('正定性: %s\n', mat2str(all(eigenvals > 0)));

% Cholesky分解
try
    L = chol(AtA_good, 'lower');
    fprintf('Cholesky分解成功\n');
    
    % 求解 L*L^T*x = A^T*b
    y = L \ Atb_good;           % 前向替换
    x_chol = L' \ y;            % 后向替换
    
    % 与直接求解比较
    x_direct_good = AtA_good \ Atb_good;
    fprintf('Cholesky解与直接解的差异: %e\n', norm(x_chol - x_direct_good));
    
catch ME
    fprintf('Cholesky分解失败: %s\n', ME.message);
end

%% 7. 迭代改进
fprintf('\n7. 迭代改进\n');
fprintf('通过迭代改进提高解的精度\n\n');

% 使用病态系统演示迭代改进
x0 = AtA_ill \ Atb_ill;  % 初始解
fprintf('初始解: x₀ = [%.6f; %.6f]\n', x0);

% 迭代改进
max_iter = 3;
x_iter = x0;

for iter = 1:max_iter
    % 计算残差
    r = Atb_ill - AtA_ill * x_iter;
    
    % 求解修正方程
    delta_x = AtA_ill \ r;
    
    % 更新解
    x_iter = x_iter + delta_x;
    
    fprintf('迭代 %d: x = [%.6f; %.6f], ||δx|| = %e\n', ...
            iter, x_iter, norm(delta_x));
end

%% 8. 加权最小二乘
fprintf('\n8. 加权最小二乘\n');
fprintf('处理不等精度观测数据\n\n');

% 权重矩阵
weights = [1; 2; 1; 3];  % 不同观测的权重
W = diag(weights);

fprintf('权重向量: [%g, %g, %g, %g]\n', weights);

% 加权正规方程: A^T*W*A*x = A^T*W*b
AtWA = A' * W * A;
AtWb = A' * W * b;

x_weighted = AtWA \ AtWb;
fprintf('加权最小二乘解: x = [%.6f; %.6f]\n', x_weighted);
fprintf('普通最小二乘解: x = [%.6f; %.6f]\n', x_normal);

% 比较残差
residual_weighted = A * x_weighted - b;
weighted_residual_norm = sqrt(residual_weighted' * W * residual_weighted);

fprintf('加权残差范数: %.6f\n', weighted_residual_norm);
fprintf('普通残差范数: %.6f\n', residual_norm);

%% 9. 性能比较
fprintf('\n9. 性能比较\n');
fprintf('正规方程法与其他方法的性能对比\n\n');

sizes = [100, 200, 500];
fprintf('矩阵规模    正规方程    QR分解    SVD分解\n');
fprintf('--------    --------    ------    -------\n');

for i = 1:length(sizes)
    m = sizes(i);
    n = min(m, 20);
    A_test = randn(m, n);
    b_test = randn(m, 1);
    
    % 正规方程法
    tic;
    AtA_test = A_test' * A_test;
    Atb_test = A_test' * b_test;
    x_normal_test = AtA_test \ Atb_test;
    time_normal = toc;
    
    % QR分解法
    tic;
    [Q_test, R_test] = qr(A_test, 0);
    x_qr_test = R_test \ (Q_test' * b_test);
    time_qr = toc;
    
    % SVD分解法
    tic;
    x_svd_test = pinv(A_test) * b_test;
    time_svd = toc;
    
    fprintf('%6d      %6.4f s    %6.4f s    %6.4f s\n', ...
            m, time_normal, time_qr, time_svd);
end

%% 10. 应用示例：多项式拟合
fprintf('\n10. 应用示例：多项式拟合\n');
fprintf('使用正规方程进行多项式拟合\n\n');

% 生成数据
x_data = linspace(0, 2*pi, 20)';
y_true = sin(x_data);
noise = 0.1 * randn(size(x_data));
y_data = y_true + noise;

% 构造Vandermonde矩阵（5次多项式）
degree = 5;
V = zeros(length(x_data), degree + 1);
for j = 0:degree
    V(:, j+1) = x_data.^j;
end

fprintf('多项式次数: %d\n', degree);
fprintf('数据点数: %d\n', length(x_data));
fprintf('Vandermonde矩阵条件数: %.2e\n', cond(V));

% 正规方程求解
VtV = V' * V;
Vty = V' * y_data;
coeffs = VtV \ Vty;

fprintf('拟合系数:\n');
for j = 0:degree
    fprintf('a_%d = %8.4f\n', j, coeffs(j+1));
end

% 计算拟合质量
y_fit = V * coeffs;
rms_error = sqrt(mean((y_data - y_fit).^2));
fprintf('均方根误差: %.6f\n', rms_error);

% 绘制结果
figure('Name', '多项式拟合 - 正规方程法');
x_plot = linspace(0, 2*pi, 100)';
V_plot = zeros(length(x_plot), degree + 1);
for j = 0:degree
    V_plot(:, j+1) = x_plot.^j;
end
y_plot = V_plot * coeffs;

plot(x_data, y_data, 'ro', 'MarkerSize', 6, 'LineWidth', 1.5);
hold on;
plot(x_plot, y_plot, 'b-', 'LineWidth', 2);
plot(x_plot, sin(x_plot), 'g--', 'LineWidth', 2);
xlabel('x');
ylabel('y');
title('多项式拟合 (正规方程法)');
legend('观测数据', '拟合曲线', '真实函数', 'Location', 'best');
grid on;

end